Bridging the size gap between density-functional and many-body perturbation theory 
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The calculation of quasi-particle spectra based on the GW approximation is extended to systems 
of hundreds of atoms, thus expanding the size range of current approaches by more than one order 
of magnitude. This is achieved through an optimal strategy, based on the use of Wannier-like 
orbitals, for evaluating the polarization propagator. Our method is validated by calculating the 
vertical ionization energies of the benzene molecule and the band structure of crystalline silicon. Its 
potentials are then demonstrated by addressing the quasi-particle spectrum of a model structure of 
vitreous silica, as well as of the tetraphenylporphyrin molecule. 
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Density-functional theory (DFT) has grown into a 
powerful tool for the numerical simulation of matter at 
the nanoscale, allowing one to study the structure and 
dynamics of realistic models of materials consisting of up 
to a few thousands atoms, these days [l[. The scope of 
standard DFT, however, is limited to those dynamical 
processes that do not involve electronic excitations. The 
most elementary such excitation is the removal/addition 
of an electron from a system originally in its ground state. 
These processes are accessible to direct /inverse photo- 
emission spectroscopies and can be described in terms of 
quasi-particle (QP) spectra 0. In insulators, the energy 
difference between the lowest-lying quasi-clcctron state 
and the highest-lying quasi-hole state is the QP band 
gap, a quantity that is severely (and to some extent er- 
ratically) underestimated by DFT Q. 

Many-body perturbation theory (MBPT), in turn, pro- 
vides a general, though unwieldy, framework for QP and 
other excitation (such as optical) spectra Q. A numeri- 
cally viable approach to QP energy levels (known as the 
GW approximation, GWA) was introduced in the 60's 
but it took two decades for a realistic application 
of it to appear and even today the numerical effort 
required by MBPT is such that its scope is limited to sys- 
tems of a few handfuls of atoms. Even so, and in spite 
of the success met by MBPT in real materials @, the 
approximations made for the most demanding of its ap- 
plications are such as to shed some legitimate doubts on 
their general applicability. This situation will be referred 
to as the size gap of MBPT calculations. 

In this letter we present a strategy to substantially re- 
duce the size gap of MBPT\ based on the adoption of 
Wannier-like orbitals 0, H, Q to represent the response 
functions whose calculation is the main size- limiting fac- 
tor of MBPT. Although we focus on QP spectra within 
the GWA, this strategy easily generalizes to optical spec- 
tra, as calculated from the Bethe-Salpeter equation Toj ] . 
Our method is benchmarked by the calculation of the 
ionization potentials of the benzene molecule and of 



the band structure of crystalline silicon, and its poten- 
tials demonstrated by case calculations on vitreous sil- 
ica and on the free-base tetraphenylporphyrin molecule 
(TPPH 2 ). 

QP energies (QPE) are eigenvalues of a Schrodinger- 
like equation (QPEq) for the so-called QP amplitudes 
(QPA), which is similar to the DFT Kohn-Sham equation 
with the exchange-correlation potential, V xc (r), replaced 
by the non-local, energy-dependent, and non-Hcrmitian 
sclf-cncrgy operator, E(r, r',E) (a tilde indicates the 
Fourier transform of a time-dependent function). Set- 
ting £(r,r';.E) = — T^rpj (p being the one-particle den- 
sity matrix) would turn the QPEq into the Hartree-Fock 
equation. The next level of approximation is the GWA 
[J] where E is the product of the one-electron propagator, 
G, and of the dynamically screened interaction, W: 

Z G w(r,r';t-t')=iG(r,r , ;t-t')W(r,r';t-t'), (1) 
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(1 — P ■ v)^ 1 ■ P is the reducible electron polarization 
propagator (polarizability), P its irreducible counterpart, 
v (r, r'; t — t') = rgzpi ^{t — t') is the bare Coulomb inter- 
action, n and V are the electron density distribution and 
external potential, respectively, and a dot indicates the 
product of two operators, such as in v ■ x(r, r', t — If) = 
f dr"dt"v(r,r";t - t") X (r" ,r';t" -t'). 

The GWA alone does not permit to solve the QPEq, 
unless G and W are known, possibly depending on the 
solution of the QPEq itself. One of the most popu- 
lar further approximations is the so called G°W° ap- 
proximation (G°W°A), where the one-electron prop- 
agator is obtained from the QPEq using a model 
real and energy-independent self-energy, such as e.g. 
E° = V xc (r)5(r — r'), and the irreducible polarizabil- 
ity is calculated in the random-phase approximation 
(RPA): G°(r,r';r) = iE t ^WCMe _ie,T 9H - 
* S c rfcirfipc (r / ) e_4ecT ^( T ) ("0 an d e are zero-th order 
QPAs and QPEs, referred to the Fermi energy, v and 
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c suffixes indicate states befow and above the Fermi en- 
ergy, respectively, and is the Hcaviside step function) 
and P°(r,r';r) = -iG°(r,r';r)G (r',r; -r). To first or- 
der in £' = T,qo W o — S°, QPEs are given by the equation: 



(2) 



where (A) n = (ip n \A\ip n ). 

The apparently simple G°W°A still involves severe dif- 
ficulties, mainly related to the calculation and manipu- 
lation of the polarizability that enters the definition of 
W°. These difficulties are often addressed using the so 
called plasmon-pole approximation [B| , which however in- 
troduces noticeable ambiguities and inaccuracies when 
applied to inhomogeneous systems A well estab- 



lished technique to address QP spectra in real materials 
without any crude approximations on response functions 
is the space-time method (STM) by Godby et al. [13]. In 
the STM the time/energy dependence of the G°W° oper- 
ators is represented on the imaginary axis, thus making 
them smooth (in the frequency domain) or exponentially 
decaying (in the time domain). The various operators 
are represented on a real-space grid, a choice which is 
straightforward, but impractical for systems larger than 
a few handfuls of inequivalent atoms. In this paper we 
combine the imaginary time/frequency approach of the 
STM with a novel representation of the response func- 
tions, based on localized Wannicr-likc orbitals, thus en- 
hancing the scope of MBPT calculations so as to embrace 
systems potentially as large as a few hundreds atoms. 

In the STM, the self-energy expectation value in Eq. 
((21) is obtained by analytically continuing to the real fre- 
quency axis the Fourier transform of the expression: 



^ n (r)Vi(r)^(r')^„(r , )W(r, r'; i T )dvdv' , (3) 



where the upper (lower) sign holds for positive (negative) 
times, the sum extends below (above) the Fermi energy, 
and QPAs are assumed to be real. By substituting v for 
W, Eq. ([3]) yields the exchange self-energy, whereas v-U-v 
yields the correlation contribution, Ecs whose evaluation 
is the main size-limiting step of GW calculations. 

Suppose that a small, time- independent, basis set 
to represent the polarizability exists: II(r, r',ir) « 
E F n^(«T)$ M (r)i„(r'). Eq. © then gives: 

(£c(*t)>„ ^TY.^^^^nhuOiEl-ei), (4) 
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r<f>„(r')drrfr' and E c is 



where S n i tV — j ym 1 j | r - r '| "^V 1 i"- LUL aw- 
an energy cutoff that limits the number of conduction 
states to be used in the calculation of the self-energy. A 
convenient representation of the polarizibility would thus 
allow QPEs to be calculated from Eq. by analytically 



continuing to the real axis the Fourier transform of Eq. 
(|4"1) . Such an optimal representation is identified in three 
steps: i) we first express the Kohn-Sham orbitals, whose 
products enter the definition of P°, in terms of localized, 
Wannier-like, orbitals; ii) we then construct a basis set 
of localized functions for the manifold spanned by prod- 
ucts of Wannier orbitals; finally, Hi) this basis is further 
restricted to the set of eigenvectors of P° , corresponding 
to eigenvalues larger than a given threshold. 

Let us start from the RPA irreducible polarizability: 



P°(r,r';z^) 



cMr)<Mr')x™M, 



(5) 



where x° cv {iuj) 



2Rc 



and $ c „(r) = 



tp c (r)ilj v (r). We express valence and conduction QPAs 
in terms of localized, Wannier-like, orbitals: 



V 

J (r)=^V-Vc(r)^(e c )^-e c ), 



(6) 



where E c < E c is a second energy cutoff that limits 
a lower conduction manifold (LCM) to be used in the 
construction of the polarization basis. According to the 
choice of the U and V matrices, the u's and v's can be 
either maximally localized 0, @] or non-orthogonal gener- 
alized [l3j Wannier functions. We then reduce the num- 
ber of product functions from the product between the 
number of valence and conduction states, which scales 
quadratically with the system size, to a number that 
scales linearly. To this end, we express the $'s as ap- 
proximate linear combinations of products of the u's 
v's: $™(r) » J2 rs O CVtrs W rs {r)6{\W rs \ 2 - 8l ), where 
O cv y v > = U VV 'V cc f, W rs (v) = u r (r)v s (r), \W CV \ is the 
L 2 norm of W rs (r), which is arbitrarily small when the 
centers of the u r and v s functions are sufficiently dis- 
tant, and si is an appropriate cutoff. The number 
of products can be further reduced on account of the 
non-orthogonality and mutual linear dependence of the 
W's. To this end, let us define the overlap matrix: 
Q pcr = J W p (r)W a (r)dr, where the p and a indices stand 
for pairs of rs indices. The magnitude of the eigenvalues 
is a measure of linear dependence, and an orthonormal 
basis can be obtained by retaining only those eigenvectors 
U v whose eigenvalue, q vi is larger than a given threshold, 
s 2 : ¥„(r) w -j^YspUvpWplT), for q v > s 2 . A final 
(and practically ultimate) refinement can be achieved by 
diagonalizing P° in the basis thus reduced and by re- 
taining only those eigenvectors whose eigenvalue is larger 
that a third threshold, S3. The result of this last proce- 
dure would depend on frequency, which would make it 
impractical. We have verified, however, that the mani- 
fold spanned by the most important eigenvectors of P° 
in the (imaginary) time domain depends very little on 
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FIG. 1: Calculated ionization potential of the benzene 
molecule (solid lines, left scale) and dimension of the polar- 
ization basis (dashed lines, right scale) versus the S2 cutoff. 
The polarization basis has been constructed with a conduction 
energy cutoff E^ = 16.7 eV (red, 100 states), E% = 28.6 eV 
(green, 300 states), and Eq = 38.3 eV (blue, 500 states). In- 
set: calculated ionization potential as a function of the overall 
conduction energy cutoff, Eq. Black line: experimental value; 
red line: fit to the calculated values (green triangles); blue line 
extrapolated value. See text for more details. 



time, which permits the use of a same basis at differ- 
ent frequencies. When the basis for P° is built out of 
orthonormal Wannicr orbitals, this last refinement does 
not result in any further improvement if S3 > S2 
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Once an optimal basis set has been thus identified, an 
explicit representation for the irreducible polarizability, 



P°(r, r, iu) = ^M<fV(r)<Mr'), 



(7) 



is obtained. By equating Eq. ((5]) to Eq. ([7|) and taking 
into account the orthonormality of the <I>'s, one obtains: 

= E T^T^xUi^iEh - e c ), (8) 



where T cv ^ — J <f> ct ,(r)<f> A1 (r)<ir. A representation for IT 
is finally obtained by simple matrix manipulations. 

Our scheme has been implemented in the Quan- 
tum ESPRESSO density functional package [TJI], for 
norm-conserving as well as ultra-soft pseudopotcn- 
tials, resulting in a new module called gww.x which 
uses a Gauss-Legendre discretization of the imaginary 
time/frequencies half-axes, and that is parallelized ac- 
cordingly. We first illustrate our scheme by considering 
an isolated benzene molecule in a periodically repeated 
cubic cell with an edge of 20 a.u. using a first conduc- 
tion energy cutoff E c = 56.7 eV, corresponding to 1000 
conduction states, and a cutoff on the norm of Wannicr 
products s\ =0.1 a.u. In Fig. 1 we display the de- 
pendence of the calculated ionization potential (IP) on 
the second conduction energy cutoff used to define the 
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TABLE I: QPEs (eV) calculated in crystalline silicon and 
compared with experimental (as quoted in Ref. Q2) and pre- 
vious theoretical results 'Thi' and 'TI12' indicate calcu- 
lations made with S2 = 0.01 and S2 = 0.001 a.u. respectively, 
while Np is the dimension of the polarization basis. 



polarization basis, E c , and on the cutoff on the eigen- 
values of the overlap matrix between Wannier products, 
S2- Convergence within 0.01 eV is achieved with a con- 
duction energy cutoff Eq smaller than 30 eV (less than 
300 states) and a polarization basis set of only ~ 400 
elements. The convergence of other QPEs is similar. 
The inset of Fig. [I] shows the convergence of the IP 
with respect to E c , which turns out to be unexpect- 
edly slow. These data can be accurately fitted by the 
simple formula IP(E C ) = IP(oo) + A/E c , resulting in a 
predicted ionization potential IP (00) = 9.1 eV, in good 
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agreement with the experimental value of 9.3 eV 
The potential of our method for large molecular system 
is illustrated by a calculation for the TPPH2 molecule 
(C44H30N4) in a periodically repeated orthorhombic su- 
percell of 75.6 x 75.6 x 26.5 a. u. Using values of 

31.1, 40.5 and 48.1 eV for E c (corresponding to 2000, 
3000 and 4000 conduction states) and E 2 C = 17.2 eV 
(corresponding to 750 conduction states), s\ = 1.5 and 
S2 = 0.01 a.u., which lead to a basis of 2797 elements, 
yields an extrapolated IP (00) of 6.0 eV, in fair agreement 
with the experimental value of 6.4 eV 18 1. 

In order to demonstrate our scheme for extended sys- 
tems [19T ]. we consider crystalline silicon treated using a 
64-atom simple cubic cell [l6| at the experimental lat- 
tice constant and sampling the corresponding Brillouin 
zone (BZ) using the T point only. This gives the same 
sampling of the electronic states as would result from 6 
points in the irreducible wedge of the BZ of the elemen- 
tary 2-atom unit cell. Our calculations were performed 
using Eq = 94.6 eV (corresponding to 3200 conduction 
states) and E c = 33.8 eV (corresponding to 800 states 
in the LCM), s± = 1.0 a.u. and two distinct values for 
S2 (0.01 and 0.001). In Tab. [I] we summarize our results 
and compare them with previous theoretical results, as 
well as with experiments. An overall convergence within 
a few tens mcV is achieved with a S2 cutoff of 0.001 a.u., 
corresponding to a polarization basis of *~ 6500 elements. 
The residual small discrepancy with respect to previous 
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FIG. 2: Electronic density of states for a model of vitreous 
silica: LDA (dashed line) and GW (solid line). A Gaussian 
broadening of 0.25 eV has been used. 



results jl_2J is likely due to our use of a supcrcell, rather 
than the more accurate k-point sampling used in pre- 
vious work. Our ability to treat large supcrcclls give 
us the possibility to deal with disordered systems that 
could hardly be addressed using conventional approaches. 
In Fig. [2] we show the QPE density of states as calcu- 
lated for a 72-atom model of vitreous silica 1|| 2(J . We 
used E c — 48.8 eV (corresponding to 1000 conduction 
states), Eq = 30.2 eV (corresponding to 500 states in 
the LCM), si = 1 a.u. and S2 = 0.1 a.u. (giving rise 
to a polarization basis of 3152 elements). We checked 
the convergence with respect to the polarization basis by 
considering S2 = 0.01 a.u. which leads to a basis of 3933 
elements. Indeed, the calculated QPEs differ in average 
by only 0.01, eV with a maximum discrepancy of 0.07 eV 
[2l| . The quasi-particle band-gap resulting from our cal- 
culations is 8.5 eV, to be compared with an experimental 
value of ~ 9 eV 22j and with a significantly lower value 
predicted by DFT in the local-density approximation (5.6 
eV). 

In conclusion, we believe that expressing density re- 
sponse functions in terms of localized basis sets will per- 
mit to extend the scope of many-body perturbation the- 
ory to large models of molecular and extended, possibly 
disordered, systems. The extension of the methodology 
presented in this paper for quasi-particle spectra to opti- 
cal spectroscopies using the Bethe-Salpeter formalism is 
straightforward and presently under way. 

We thank C. Cavazzoni for his help in the paralleliza- 
tion of the code. This work has been partially funded 
under the Italian CNR-INFM Seed Projects scheme. 
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